How the Interplay among Conformational Disorder, Solvation, Local, and Charge-Transfer Excitations Affects the Absorption Spectrum and Photoinduced Dynamics of Perylene Diimide Dimers: A Molecular Dynamics/Quantum Vibronic Approach

In this contribution we present a mixed quantum-classical dynamical approach for the computation of vibronic absorption spectra of molecular aggregates and their nonadiabatic dynamics, taking into account the coupling between local excitations (LE) and charge-transfer (CT) states. The approach is based on an adiabatic (Ad) separation between the soft degrees of freedom (DoFs) of the system and the stiff vibrations, which are described by the quantum dynamics (QD) of wave packets (WPs) moving on the coupled potential energy surfaces (PESs) of the LE and CT states. These PESs are described with a linear vibronic coupling (LVC) Hamiltonian, parameterized by an overlap-based diabatization on the grounds of time-dependent density functional theory computations. The WPs time evolution is computed with the multiconfiguration time-dependent Hartree method, using effective modes defined through a hierarchical representation of the LVC Hamiltonian. The soft DoFs are sampled with classical molecular dynamics (MD), and the coupling between the slow and fast DoFs is included by recomputing the key parameters of the LVC Hamiltonians, specifically for each MD configuration. This method, named Ad-MD|gLVC, is applied to a perylene diimide (PDI) dimer in acetonitrile and water solutions, and it is shown to accurately reproduce the change in the vibronic features of the absorption spectrum upon aggregation. Moreover, the microscopic insight offered by the MD trajectories allows for a detailed understanding of the role played by the fluctuation of the aggregate structure on the shape of the vibronic spectrum and on the population of LE and CT states. The nonadiabatic QD predicts an extremely fast (∼50 fs) energy transfer between the two LEs. CT states have only a moderate effect on the absorption spectrum, despite the fact that after photoexcitation they are shown to acquire a fast and non-negligible population, highlighting their relevance in dictating the charge separation and transport in PDI-based optical devices.


S1.1 Definition of intermolecular dimer coordinates
The intermolecular dimer coordinates employed in this work, and shown in Figure 3.a in the main text, are based on three unit vectors defined with respect of each PDI aromatic core. In fact, as displayed in Figure   S1, for each core k,û || lies parallel to the PDI core's long axis,û ⊥ is the unit vector normal to the aromatic plane andû S is obtained by the vector product between the first two. According to these molecular reference " PDI-1 PDI-2 Figure S1: Definition of unit vectors on each monomer. systems, the intermolecular coordinates shown in Figure 3 can be defined as follows: Figure S2: A cartoon representing the the four coupled diabatic states defined for the PDI dimer in terms of the individual monomers. The two green lines indicate the HOMO (lower line) and LUMO (upper line) of each monomer, and the arrows represent the electrons with their spin.

S1.3 Spectral intensity and first moment
Analytical expressions for the total intensity (moment 0, M 0 ) can be straightforwardly obtained considering the absorption lineshape, (ω)/ω neglecting the prefactor C , and integrating over the frequency.
where µ g,L is the transition dipole of the diabatic states representing the local excitations. The first moment, i.e. the center of gravity of the spectrum is

S1.4 Constrained dimer optimization
According to the original Ad-MD|gVH scheme, S1 the reference geometry in the reduced space of the fast coordinates, q r o , is computed performing an harmonic extrapolation from the gradient and Hessian matrix.
Here, since in the simplified route of the Ad-MD|LVC scheme specific normal modes for each snapshot are not needed, the time consuming computation of the Hessian matrices is not longer required, and q r o can be retrieved simply by re-optimizing the PDI core of each monomer. However q r,α o is still needed to perform the single-point FrD-Ex diabatization. Concretely, to retrieve q r,α o for each MD snapshot α, we carried out a DFT geometry re-optimization for each individual monomer, where the atoms pertaining to its lateral pendants and to the other PDI monomer were constrained, hence only the atoms of one PDI core were allowed to relax. An ONIOM QM/MM partition was followed, where the surrounding chemical species located within a radius of 4Å (ACN, Cl − , and the other PDI) are included in the DFT optimization at MM S4 level (employing QMD-FF point charges and LJ terms), whereas all the chemical species within a radius of 20Å are instead treated as QMD-FF point charges. Although approximated, this is the easiest procedure to perform the required constrained optimization in the space of the stiff coordinates only. Finally, for each snapshot, the optimized geometry of the dimer is straightforward obtained by combining the two optimized monomers.
In principle, to have a more accurate structure, the aforementioned re-optimization could be carried out considering both PDI monomers at the QM level, hence including the effect of the mutual polarization.
Yet, to avoid to alter the soft intermolecular coordinates describing the relative disposition of the PDI cores sampled by MD, several constraints between the PDI cores should be applied, with a significant impact on the computational cost. In order to verify the small impact of neglecting such mutual polarization, we compared the structures of the PDI dimer for 4 snapshots (2 extracted from the MD of the anti isomers and 2 from the one of the syn isomers) obtained by including the second PDI monomer at either MM or QM level. The differences between the two structures were found extremely limited, being the averaged Root Mean Square Deviation (RMSD) between them of the order of 0.06Å. To further corroborate the validity of the approximation, we performed TD-DFT calculations and adiabatic-to-diabatic transformation considering both the constraint optimization geometries and the fully-relaxed ones, to eventually compute the distance between the diabatic Hamiltonian matrices utilizing the Frobenius distance: The distances thus calculated have been averaged for the 4 snapshots considered in this analysis obtaining very small both distance and standard deviation: 0.063 ± 0.016 eV, and hence validating our procedure.

S1.5 MCTDH settings
Nuclear wavepackets were propagated for 100 fs in steps of 0.5 fs. For the parametric study and the dimer minimum calculation, we used a combination of a 7th order Bulirsch-Stoer integrator with a stepsize of 2.5·10 −4 and a 5th order SIL integrator, both with an accuracy of 10 5 . The mean field matrices were assumed constant during the time interval (Constant Mean Field, CMF), also with an accuracy of 10 5 . We used effective modes generated with a hierarchical representation of the Hamiltonian. For the convergence of the spectrum broadened with a Gaussian with HWHM=0.05 eV 2 blocks were necessary (8 coordinates when only λ ii coupling are included, 12 for λ ii and λ ij ). For the calculation of the spectra of each snapshot we adopted a narrower broadening with HWHM=0.03 eV. Therefore the effective modes of one additional hierarchy blocks were required, for a total of 12 modes. In these cases we used the Multi-Layer extension S5 ML-MCTDH with a 7th order Runge-Kutta integration and the mean field matrices were updated at each integration step (Variable Mean Field, VMF). The same ML-MCTDH settings were used for the computation of the time-evolution of the populations. Figure S3: Combination of degrees of freedom for the MCTDH calculations of the full LVC Hamiltonian, including all λ ii and the linear couplings λ 12 and λ 34 . All E ij (0) terms were included. We adopted 2 hierarchy blocks for a total number of 12 coordinates. Figure S4: ML-MCTDH tree used of the simplified route employed for the snapshot calculations, where the LVC Hamiltonian only includes the diagonal gradients λ ii and all the constant terms E ij (0). We adopted 4 hierarchy blocks for a total number of 16 coordinates for testing, 3 blocks for narrow spectra (HWHM = 0.03eV) and 2 hierarchy blocks for a total number of 8 coordinates for broader spectra (HWHM=0.05eV). The tree for 2 hierarchy blocks was obtained by suppressing the right branch of coordinates, while on the tree for 3 blocks, the right branch of coordinates was suppressed the coordinates were combined in groups of 3 on the lowermost level. Figure S5: Convergence test for the number of hierarchical blocks. Spectra obtained using the LVC gradients λ ii and the constant coupling terms E ij (0) calculated for the minimum of the PDI dimer. The spectra were redshifted by 0.25eV and convoluted with a Gaussian with HWHM = 0.03eV. The calculations with 3 and 4 blocks include respectively 12 and 16 coordinates.

S2.1 Experimental spectra
Difference spectra were obtained by subtracting, from the initial spectrum of PDI in pure acetonitrile, the subsequent spectra in the presence of increasing water concentration ( Figure S7) or the spectrum obtained by adding 20µl of the saturated solution on PDI in ACN to water solution ( Figure S9). In order to do so, each spectrum was normalized between 0 and 1 by using the lowest wavelength absorption peak and by subtracting the background baseline between 600 and 800 nm. In fact, the spectral baseline progressively increased upon the addition of water due to some light scattering caused by aggregate formation. Difference spectra (ACN) 0.5ml 1ml 1.5ml 2ml Figure S7: Difference between the spectra registered for PDI(ACN) in the presence of different amount of water (0.5 ml -2 ml) and the spectrum registered for the saturated solution of PDI in ACN (black line in Figure S6).  Figure S9: Difference between the spectrum registered when 20µl of the PDI(ACN) solution is added in water and the spectrum registered for the saturated solution of PDI in ACN (black line in Figure S6).
As discussed in the main text, the full LVC protocol was employed to perform a static computation of the vibronic spectrum in ACN solution. With the term "static" we mean that we neglected all the effects of the fluctuations of the soft coordinates and we computed the spectrum considering a single reference geometry: the equilibrium geometry of the dimer in the anti conformation, obtained by a full QM optimization in ACN (C-PCM S2 ). The coordinates of the flexible lateral chains, i.e. all bonds, angles and dihedrals related to atoms in the chain, which are weakly involved in the electronic transition, were projected out of the coordinate space, following the iterative procedure introduced in previous works. S1,S3 The inter-monomer coordinates, defined adopting the translation-rotation-internal coordinate (TRIC) system introduced by Wang and Song, S4 were also projected out using the same protocol. This reduces the coordinate space from 450 normal coordinates to 224. The LVC parameters of the reduced dimensionality system were thereafter computed as described in the Computational Details. This step is computationally costly as it requires 449 TD-DFT computations of the dimer (152 atoms).
Thereafter, the absorption spectrum was eventually computed adopting a hierarchical representation by blocks of the Hamiltonian. S5-S9 Two different sets of collective coordinates were tested for the dimer minimum. The first set included all the relevant LVC parameters (gradients and linear coupling) while the second contained only the gradients. Figure S3 (2 blocks, 12 collective coordinates) and Figure S4 (2/4 blocks, 8/16 collective coordinates) give details on the adopted ML-MCTDH tree. We also tested the convergence of the spectra with respect to the number of hierarchical blocks in Figure S5. Both spectra are very similar and for larger broadening they are identical, however, for HWHM=0.03eV 3 blocks do not fully converge the high energy region and 4 blocks most be used.

S2.3 Accuracy of the simplified Ad-MD|gLVC route
In this section we test the accuracy of the approximations that allow us to speed up the full computational protocol presented in the left panel of Figure 2, leading to the simplified approach shown in the right panel of the same Figure. To this end, Figure S10 shows the results of a static computation of the vibronic spectrum in ACN solution. The computed spectrum is reported in black in Figure S10, and is only marginally different from a spectrum computed with λ ij = 0 (green line), proving that this approximation does not introduce any significant error. This optimal result is not surprising considering that the constant coupling terms E ij are much larger than the norm of the vector of the corresponding linear terms |λ ij |, as evident from Table A.
It is also noteworthy that, with this simplification, the number of coordinates corresponding to two blocks of the hierarchy is smaller (8) since it is not necessary to include coordinates parallel to the λ ij , making the MCTDH propagations faster. This is especially relevant for those snapshots where calculations with 4 blocks are required for high-resolution spectra: omitting λ ij reduces the dimensions of the vibrational space from 24 to 16.
We also tested the second approximation proposed in the right panel of Figure 2, namely that the diagonal gradients λ ii may be computed on monomers, rather than the dimer. Despite this crude approximation,  S3 The only appreciable difference is that the vibronic peaks of the approximated spectrum appear slightly more resolved, but this effect is then washed out by the average over a representative set of configurations of the dimer.

S2.4 MD Structural Analysis: lateral pendants
Figures S11-S13 show the distribution and the time evolution of the dihedral angle defining the rotation of the lateral pendants with respect to the π-plane of the PDI core (δ1). To highlight the rotation around δ1, and therefore the inter-conversion from the anti/anti dimer to the syn/anti one, we defined the starting value of this dihedral angle as always positive. Such definition allows us to easily recognize the anti → syn transition because, when this occurs, δ1 assumes negative values.

S14
S2.5 MD Structural Analysis: aggregate shape in ACN Figure S14 shows the time evolution of the distance between the center of the PDI cores (see ρ in Equation S1 and in Figure 3 of the main text). This Figure shows the tendency of PDI to form self-assembled dimers and has been computed along the MD simulations of both the isomers (syn and anti) for the PDI 2 @ACN system. The 2D colour plots (Figures S14-S20) display the cross-correlation of the distribution of the intermolecular geometrical parameters analyzed in this work. The color scale reported on the right hand side of these graphs represents the frequency with which it is populated a PDI 2 structure having the intermolecular geometrical characteristics displayed in the X and Y axes. Also in this case, as for Figure S14, the MD simulations have been analyzed of both isomers of the PDI 2 @ACN system.        Figure 3 and Figure S1. S18 S2.6 MD Structural Analysis: aggregate shape in H 2 O Figure S21 shows the distributions of the intermolecular descriptors used to analyze the dimer shape in ACN reported in Figure 4 of the main text. Recalling that the anti isomers only have been used for the MD simulations in water and that the anti → syn transition is faster than in ACN (see Figures S11-S13), the distributions of Figure S21 do not differ substantially from those computed for the corresponding simulations in ACN (Figure 4a). The only difference worth mentioning concerns the β angle distribution that shows smaller values in Figure S21, demonstrating the tendency of PDI to form dimers with more parallel planes in aqueous solutions.  Figure S21: Distributions of inter-molecular (PDI-PDI) angles, distances, and displacements as defined in Figure 3 and Figure  S1 computed along the MD trajectories of the PDI 2 @H 2 O starting with the anti isomer. S19 S2.7 Parametric study of stacked dimers Figure S22 complements Figure 5 in the main text. While for the sake of brevity, in the latter figure we only reported the energy gap between the CT and LE states, here we plot the absolute energies of the four diabatic states. It is interesting to notice that LE energies change only very moderately. On the contrary CT energies decrease remarkably with the three distances and in particular with R π and then R S . They are insensitive to α and more sensitive to the two β angles, which also induce a difference in energy between the two CT states. S20 Figure S23 shows the absorption spectra computed for the dimer configurations investigated in the parametric study. Notice that spectra were red-shifted by the same amount (0.27 eV) adopted to plot the spectra in Figure 6 in the main text. It is evident that R decreases upon decreasing either one of the three considered distances, or α. In all these cases, these movements correspond to an increase of the exciton coupling. On the contrary, the spectral shape seems rather insensitive to β(R S ), whereas large tilting around β(R ) lead to a number of weak bands on the red wing of the spectrum. Figure S23: Spectra of the PDI dimer without pendants and in gas phase for different structural parameters. All calculated spectra were redshifted by 0.27eV and arbitrarily scaled so to fit in the same figures. We applied the same scaling factor to all spectra in the same panel, but such factor can vary from panel to panel. Spectra convoluted with a Gaussian of HWHM = 0.05 eV.

S2.8 The prediction of different models for PDI monomer spectrum
It should be noticed that the monomer spectrum reported in Figure 6 in the main text is slightly different from the Ad-MD|gVH one given in Ref. S3. Indeed, to be coherent with the approximations in the LVC model, here we neglected quadratic differences in the PESs from the ground to the excited states, i.e.
Duschinsky rotation and changes in the vibrational frequencies, replacing the VH model by the simpler VG one. S10 Ad-MD|gVH and Ad-MD|gVG are compared in Figure S24, showing that such quadratic corrections (neglected also in the spectra of the dimer) have a very slight impact on the spectral shape and simply cause a red-shift by ∼ 0.08 eV. Figure S24: Comparison of experimental and calculated spectra for the PDI monomer. The calculated spectra were obtained either from the dimer snapshots, setting all inter-state couplings to zero (PDI 2 @ACN Ad-Md|gVG) or by applying the Ad-MD|gVH protocol whose data are reproduced from a previous publication, S3 or, finally, still adopting the same electronic data of ref. S3 but employing the simpler Ad-MD|gVG protocol. The latter model, like the LVC one, assumes that all excited electronic states have the same normal modes and frequencies of the ground state. All calculated spectra were redshifted by 0.27eV and arbitrarily scaled by the same factor. The stick transitions were convoluted with a Gaussian of HWHM=0.05eV. Figure S24 also compares the monomer spectrum with Ad-MD|VG and with Ad-MD|gLVC, setting to zero all inter-state couplings (PDI 2 @ACN Ad-MD|gVG). This is equivalent to compute the spectrum of one monomer in a solvation shell comprising also the second monomer. The position of the spectrum is practically S22 the same, showing that the effect of one PDI on the vertical excitation of the other monomer is very small. The intensity of the second peak is moderately larger for the Ad-MDg|LVC, and this finding is attributed to the effect of the replacement of the pendants by H atoms (notice that a similar effect was appreciated in Figure S10). Comparison with the experiment shows that all the computed spectra overestimates the ratio R and the larger discrepancy is observed for the "PDI 2 @ACN Ad-MD|gVG" calculation. As discussed in the main text, this suggests that the simplified model adopted to compute the gradients of the LE states can contribute to the excessive splitting of the lowest two bands in the dimer spectrum.
S2.9 Further investigations on the origin of the lowest energy peak underesti-

mation. A small two-states two-modes model
To gain insight on the origin of the spectral bands, we performed time-independent calculations of the spectra on a fully eclipsed PDI dimer (R π = 4.0Å). In this framework, including the nonadiabatic effects is very computational demanding so we used a very reduced model for the dimer composed of 2 electronic states (the two LE) and two normal coordinates (the one carrying the largest gradient for each state). In addition, the constant interstate coupling E 12 (0) was also included in the model while the linear vibronic coupling terms λ ij (k) were neglected, like it is done in the simplified route adopted in the paper. Moreover the coupling to the CT states is neglected. The analysis of Figure 7 shows that it has a very limited impact on the spectral shape. Whereas the two states are the two local excitations |L 1 and |L 2 , the two considered coordinates are the first members of the hierarchy, i.e. the coordinates q 1 and q 2 describing the displacement of the L 1 and L 2 equilibrium geometries with respect to the ground state geometry.
First, we check convergence by comparing with the spectrum obtained for the same model in the timedependent framework for which convergence is easily achieved. Figure S25 shows that the TI spectrum agrees very well with the TD one so we can confidently assign the spectral bands from this model. Figure S26 compares the predictions of the two-states two-coordinates model (2Q) with those of the model adopted in the main text and involving all the normal coordinates (All Q). We report both nonadiabatic spectra (LVC) of the dimer and the diabatic limit obtained switching off the excitonic coupling and corresponding to the monomer spectrum (FC|VG). Focusing first on the monomer spectra, the figure shows that the "2Q" model reproduces the general shape of the spectrum obtained including all coordinates, but for a relevant overestimation of the ratio R. Switching on the excitonic coupling, the 2Q model nicely captures the change of the spectrum from the monomer to the dimer predicted by the "All Q" model, thus reproducing the main experimental observations. Like for the monomer case, also for LVC nonadiabatic spectra of the dimer the "2Q" model predicts larger R than the "All Q" model. Figure S25: Comparison of the calculated LVC time-independent and time-dependent spectra for a reduced dimensionality model including 2 normal coordinates and 2 excited states. Stacking distance of 4Å and relative rotation angle α = 0°. For TD approach, both auto and cross terms were included in the computation of the spectrum. Figure S26: Calculated LVC spectra for the two-states-two-coordinates simplified model and the all coordinates model turning on or off the interstate coupling between the local excitation. The CT states were excluded of the model. To simplify the comparison, all spectra have unit area. Stacking distance of 4Å and relative rotation angle α = 0°.

S24
After proving that the 2Q model capture the main features of the monomer and dimer spectra, we employ it to further analyse the role of the exciton coupling and the gradient of the LE state (which determines the displacement of its equilibrium). For such an analysis the 2Q model is quite helpful because it is so small to allow us the computation of the vibronic eigenstates whose assignment is straightforward introducing symmetry-adapted electronic states |A 1 = 1/ √ 2(|L 1 − |L 2 ), |A 2 = 1/ √ 2(|L 1 + |L 2 ) and coordinates q − = 1/ √ 2(q 1 − q 2 ), q + = 1/ √ 2(q 1 + q 2 ). With reference to a perfectly cofacial arrangement, the symmetric |A 2 is higher in energy and bright, whereas the anti-symmetric one |A 1 is lower in energy and dark. Figure S27: Calculated LVC time-independent spectra and assignment of the vibronic bands on the diabatic representation (top) and the symmetry-adapted adiabatic representation (bottom) for a two-states-two-coordinates simplified model. In the diabatic representation each label (1 and 2) correspond to one monomer, while for the symmetry-adapted adiabatic representation the states are ordered in increasing energy and since the excitation coupling is positive, S1 and S2 correspond to the antisymmetric and symemtric combination of the local excitations. The coordinates q + and q − are respectively the symmetric and antisymmetric combinations of the monomer normal coordinates. "0" stands for the 0-0 vibronic transition. Stacking distance of 4Å and relative rotation angle α = 0°. For this arrangement S1 is totally dark and S2 is bright. Figure S27 shows the assignment of the three bands in terms of the normal coordinates on the monomers (i.e., q 1 and q 2 , top panel) and also in terms of symmetry adapted coordinates q + = q 1 + q 2 and q − = q 1 − q 2 S25 (bottom). Before starting the analysis, it is worth stressing that the position of the bands of the 2Q spectrum does not correspond to the one of the spectra reported in Figure 6. This finding is due to a number of reasons (i) the fact that the spectra in Figure 6 were red-shifted by 0.27 eV, (ii) the different spatial arrangement of the two monomers (an average from MD snapshots in Figure 6, a cofacial arrangement with a distance 4Å here), (iii) the lack of the effect of the other coordinates. Therefore the 2Q spectrum has three bands at ∼ 2.78, 2.95 and 3.12 eV, which correspond to the three bands in Figure 6 at ∼ 2.3, 2.5 and 2.7 eV.
The lowest two bands of the spectrum (those at ∼ 2.78 and ∼ 2.95 eV, respectively) are combinations of the two states |A 2 , 0 − , 0 + (bright) and |A 1 , 1 − , 0 + (dark), where the first label specifies the electronic state and the second and third numbers give the quantum number of the vibrations q − and q + . The highest band at ∼ 3.12 eV arises from the maximum band with a progression along the total-symmetric mode q + and is a combination of |A 2 , 0 − , 1 + and |A 1 , 1 − , 1 + . All the data for the assignment of the bands are given in Table B. Table B: Contribution of vibronic states to the first (top), second (middle) and third (bottom) bands upon variation of E 12 (0) (third column, λ 12 (q) fixed at -0.13 eV) and the linear coupling term λ 12 (q) (last column with E 12 (0) fixed at 0.126 eV). The states are labelled as (e, nq − ,nq + ) where e represents the adiabatic electronic state and nq − ,nq + respectively are the quanta on the coordinates q − ,q + . Finally, a parametric study changing the exciton coupling E L1,L2 and the gradient λ along the two coordinates q 1 and q 2 displayed in Figure S28 shows that the ratio of the intensity of the central band at 2.95 eV with respect to the lowest-one at 2.78 eV increases at the increase of E L1,L2 . At variance, it decreases at the increase of λ, but this is an indirect effect due to the transfer of intensity to the third band at 3.12 eV.
The splitting of the two lowest bands increases both with E L1,L2 and λ. This suggests that the disagreement in the relative intensity of these bands with respect to the experiment should be due to an overestimation of the exciton coupling. This latter contributed also to the excessive splitting of the two bands, in this case in combination with the overestimation of the displacement of the equilibrium positions of the local excitations (already noticed for the monomer). Figure S28: Calculated LVC time-independent spectra for a two-states-two-coordinates reduced dimensionality system with different values of the local excitation interstate constant coupling E 12 (0) (top, λ 12 (q − ) fixed to -0.13eV) or the linear coupling λ 12 (q) (bottom E 12 (0) fixed to 0.126eV). Therefore, the blue spectrum is the same in both panels. All stick transitions were convoluted with a Gaussian of HWHM = 0.05eV and computed spectra were redshifted by 0.45eV. The experimental spectrum is added as reference. Stacking distance of 4Å and relative rotation angle α=0°.

S27
S2.10 Comparison of the Ad-MD|LVC spectra for the PDI dimer in for the individual snapshots in ACN and water Figure S29: Ad-MD|LVC averaged spectrum and the individual contribution from each snapshot for PDI 2 anti conformer at ACN (top) and water (bottom). All stick transitions were convoluted with a Gaussian of HWHM = 0.03 eV. Table C: Interstate coupling, in eV, obtained with different functionals, all in combination with the 6-31G(d) basis set. PDI dimer with stacking distance of 4Å and relative rotational angle α = 0°. For comparison the excitonic coupling between the two local excitations has been computed also from the Coloumb interaction betweeen transition densities (EET method in Gaussian 16, using the CAM-B3LYP functional S11 ). Labels for excited states are 1 = |L 1 , 2 = |L 2 , 3 = CT (1 → 2) and 4 = CT (2 → 1)

S2.12 Relation between aggregate structure and excited states couplings
To have an accurate description of how the PDI dimer shape influences the exciton nature and delocalization, for all snapshots extracted from MD simulations in ACN and used to compute the spectra, we analyzed the Diabatic Hamiltonian matrices and displayed the relations between the geometrical descriptors used throughout the manuscript and the interstate couplings (off-diagonal terms of diabatic Hamiltonians) through 2D color maps. Whilst the color scales in Figures S15-S20 represent the probability of finding the PDI dimer in a specific geometrical configuration, the color scales in Figures S30-S41 represent the magnitude in eV of the interstate coupling between local excitations (LE-LE Coupling) and between the |L 1 and |CT (1 → 2) states (LE-CT Coupling). Although in the diabatic Hamiltonian there are four possible LE-CT couplings terms, these do not differ substantially in magnitude, therefore for the sake of simplicity, the coupling between |L 1 and |CT (1 → 2) only has been considered.            at position (0) taken over the snapshots extracted from anti and syn trajectories in acetonitrile. It clearly shows that for syn trajectory vertical energies are slightly smaller (with smaller standard deviation for CT states) and couplings slightly larger (with larger standard deviations). Table E reports the values of the elements of the LVC Hamiltonian matrix at the reference position (0) for a specific snapshot in acetonitrile where the CT(1→2) and CT(2→1) populations? are respectively very large and vanishingly small. By repeating the computation of the matrix after removing all solvent molecules we clearly show the remarkable effect played by the specific solvent configuration. In gas phase in fact the two CT states are almost degenerate and less stable than the local excitations by ∼ 0.2 eV. On the contrary, in solution CT(1→2) and CT(2→1) are respectively stabilized and de-stabilized by ∼ 0.5 eV, so that the former becomes strongly more stable than the initial state of the dynamics L 1 while the latter is so higher in energy that cannot be populated. S2.14 Absorption spectrum and population dynamics predicted by the averaged Hamiltonian Figure S42: Calculated spectra obtained by averaging the individual spectra of each snapshot (red) and the spectra obtained from a wavepacket propagation under the effect of a single LVC Hamiltonian with values obtained by the average of E d ii (0) and E d ij (0) values over the set of 100 snapshots (green and blue). Notice that we also added the results with a larger broadening (blue) selected so try to phenomenologically account for the effect of the fluctuation of the parameters, which is lost in the average procedure. Average values taken from Table D, anti conformer PDI dimer in ACN. Figure S43: Population dynamics of the diabatic states after photoexcitation to |L 1 obtained by averaging the populations of 100 snapshots (red) or from a single propagation of a single LVC Hamiltonian with values obtained by the average of E d ii (0) and E d ij (0) values over the set of 100 snapshots (green). Average LVC values taken from Table D, anti conformer PDI dimer in ACN.